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This paper studies the applicability of the 
path integral solution technique for esti- 
mating extreme response of nonlinear dy- 
namic oscillatOTS whose equations of mo- 
tion can be modelled by the use of Ito 
-Stochastic differential equations. The state 
vector process associated with such a 
model is generally a diffusion process, 
and the prtAabilily density function of the 
state lector thus saii^fie^ ihe Fokker- 
Planck-Kolmoforov equation. U is shown 
that the path integral solution technique 
combined with an appropriate numerical 
scheme constitutes a powerful method for 
solving the Fokkcr-Planck Kolmogorw 
equation with natural boundary condi- 



tions. With the calculated probability den- 
sity function of the state vector in hand, 
one can proceed lo calculate the required 
quaniities for estimating extreme re- 
spfflise. The proposed method distinguishes 
itself by remarkably high accuracy and 
numerical robustness. These features are 
highlighted by application to example 
studies of nonlinear oscillators excited by 
white noise. 
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1. Introduction 

An important element in the safety assessment of 
many engineering systems, is the task of estimating the 
probability of extreme events that may jeopardize the 
structure in some specified sense. Very often, this prob- 
lem can be formulated as finding the probability that 
some time varying random quantity does not exceed a 
specified capacity level during a given time period. 
Stated this way, the problem typically reduces to a study 
of the extreme values of a stochastic process originating 
as the response of a system subjected to some stochastic 
loading process. 

In this paper the focus will be on the problem of 
estimating the extreme response of nonlinear dynamic 
systems subjected to random forcing processes. In re- 
cent years the methods of time domain Monte Carlo 
simulations, see, e.g., Refs, [1-51, have received consid- 



erable attention as a tool for estimating response statis- 
tics. These methods are versatile and attractive in the 
sense that nonlineariiies can be easily dealt with. The 
main drawback at present is the large CPU times needed 
for accurate prcdictbn of extreme responses. Even if 
this issue seems to become less of an obstacle every 
year, portending perhaps that such methods may domi- 
nate practical estimation of response statistics of nonlin- 
ear systcm-s in the not too distant future, it will .still be 
desirable to h<i\'s available alternative methods of calcu- 
lating the response statistics, both simplified and more 
elaborate. Here we shall explore a method ba.sed on the 
theory of Markov diffusion processes. The justification 
for using this theory is related to the fact that the re- 
sponse of nonlinear dynamic systems to broad band 
random excitation can very often be accurately dc- 
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scribed by applying the theory of multidimensional 
Markov processes. By this, the extensive theory of 
Markov diffusion processes can be brought to bear on 
these problems. In particular, it can be shown that the 
probability law of response quantities can be derived by 
solving a partial differential equation, viz.. the Fokker- 
Planck (-Kolmogorov) (FPK) equation, see Refs. [6,7]. 
In most cases of practical interest, this equation has to be 
solved nunierically. 

in the next section we shall describe a method for 
solving the FPK equation that is based on a formal 
solution of the same equation. This solution is obtained 
by invoking the fact that a Markov diffusion process 
locally looks like a Browman motion. By using the 
Markov property, the global solution can then be con- 
structed by linking the local solutions, which are known 
explicitly. The obtained solution is generally known as a 
path integral solution (PIS). The reader is referred to 
Ref. [7] for a further discussion. One of the first efforts 
to exploit the PIS method explicitly in developing nu- 
merical solution algorithms is described in Ref. {$]. 
Subsequently, other authors have also used the PIS ap- 
proach to solve various random vibration problems, cf. 
Refs. [9-13]. 

Before embarking on a description of the PIS method, 
it is expedient to briefly show h(Tw the obtained solutions 
arc used in an extreme value analysis. Assuming that the 
response quantity of interest is a scalar (real) stationary 
stochastic process, Z{t) say, the PIS method typically 
provides a numerical estimate of the joint probability 
den.sity function (PDF)/zz(*,*) of Z(f) and Z(f)=dZ(r)/ 
df . It is now assumed that the mean level upcrossing rate 
Vz(,') of Z(f) can be calculated from Rice's formula as 
folk>ws 






ii{zy-\yfzz(z,y)(iy. 



(1) 



Adopting the assumption that upcrossing of high levels 
are statistically independent events, which leads lo Pois- 
son distributed crossings, it follows that an asymptotic 
approximation of the probability distribution function of 
the extreme value of the process Z(/) during a time T, 
denoted by M{T) (-sup{Z(r):Osf<r}). is given by 

Prob{M{T)^z}-cxp{-i$Xz)T} (T^^). (2) 

The accuracy of Eq. (2) depends to a large extent on 
the effective bandwidth of the response process Z(f). 
Decreasing bandwidth leads c\'cntually to a significant 
clumping effect of targe response peaks, invalidating 
the assumption of statistically independent upcrossing of 
high levels. Methods that aim at correcting for this effect 
have been proposed for Gaussian (Refs. [14,15]) and 



non-Gaussian (Ref. [16]) processes. However, this point 
will not be pursued any further here. We shall assume 
that Eq. (2) provides an acceptable approximatkjn. 
Hence, the central parameter to be determined is the 
upcrossing frequency i^K*). which Ls easily calculated 
once the joint PDF/zz(*,') of Z(0 and Z(r)-dZ(r)/dr has 
been made available. In the next section it is shown how 
this PDF can be calculated for the response of a wide 
range of nonlinear oscillators subjected to white noise 
or filtered white noise loading. 



2. The Path Integral Solution 

The path integral solution (PIS) method is suitable for 
calculating the joint probability density function (PDF) 
of a vector process X{t)=lX\{t)....Ji„{t)f (T-transposi- 
tion) satisfying a stochastic differential equatkm of the 
folkiwing form, cf. Ref. [6]. 



dX(.i)-m[X{t)]dt+Q[X{t)\dW(l). 



(3) 



Here m(*)»[m((-),...,ffi„(')]^ ntji*} denotes a real func- 
tion of n real variables. Q{')-{qtji')) denotes an nXm- 
inatrix where each q,j{') is a real function of n real 

variables. W(f)=[W,(r) W„(t)V where Wj{t), 

j^ I ,.../« are standard, real Brownian motkin processes, 
which are mutually independent, see e.g., Refs. [6,7], 
That is, £[Wj(/)J=0 and 



fc[dVV,(r)dlVj((+T)]-5„8,,,.,dT, (^l,...,w. 



(4) 



where &^.-l for x^y. 8j,=0 for xj^y. Equation (4) 
is a short-hand notation for ihe relation 

E[\\h{s,i)d\Vi[s)dWj{t}]-^Jh(,t,t)(it, where h(:') 

is a non-random function. 

Equatron (3) is interpreted here a.s an Ito stochastic 
differential equation (SDB). Since it is often relevant to 
consider Eq. (3) as being obtained as a limit of equations 
with band limited noise processes, it may happen that 
m{*) should contain correction terms to ensure a consis- 
tent limiting solution, cf. Ref. [6]. It is assumed here that 
this consideratron has already been made, and that Eq. 
(3) has the final form to be used subsequently. 

It is demonstrated in Ref. [6] that the solutkin Xit) to 
Eq. (3) is a Markov vector process. Its transitk>n proba- 
bility density function (TPD), p(x.r I j'.r'), is defined 
by the equation 

Prob{X{t)GA\XU')'X'}=\"-\p{x,t\x',t')ilx, (5) 

where /i C R" is some event, jc, x' G R", djc=dvi...dr„. 
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Provided that m{') and Q(') satisfy certain regularity 
conditions, see Rcf. [6], it can be proved that the TPD 
p(x,t \x\i') (f^r'aO) is the solution of a partial differ- 
entia! equation of the form 



O denote appropriate zero-matrices and G(*) denotes an 
(n— r)X(rt— r)-matrix function with elements g,>(*), 
/j-r+I,....M. G(*) will be called the reduced diffusi<in 
matrix. Equation (6) can now be rewritten as 



|p(x,r 1jc'/)=-X^K(J«^)p(^.' U'.'')] 



1^^ B' 



+ \l.I.r~^\g,ix)pix,t\x\n], (6) 

where G{x)^(^,{x))=Qix)Q{xf=il.l^<l^qA ^ 

with initial condition p(x,/' |jr',f')=5(x-j:'). G(-) will 
be called the diffusion matrix and Bq. (6) will be re- 
ferred to as the Fokker-Planck-Kolmoeorov (FPK) 
equation. Since clearly Prob{X(t)&R''\X{i')^x'}~\ 
the TPD satisfies the following normalization condition 



(••' /?(x,/U',r')dj:=l. 



(7) 



Let/(jc,/) denote the PDF of the random vector A'(r). 
If f(x,i')=w{x) for some initial PDF w(x), then it is 
recognized from Eq. (6) and the relation 



/(x,r)=J— Jp(r,r U'.r'>f (x')djr' 



(8) 



lhat/(x,/) itself is a solution of Eq. (6) satisfying the 
initial cond ition /( J ,/')-"♦■ (x). 

In this paper wc shall be interested primiuily in sta- 
tionary solutions /,(x) to Eq. (6), that is 



/,(x)=lim/(x./)=!imp(x,r !x',f') 



(9) 



provided they exist. Even when both limits exist, it is 
clear that lim/(x,/) provides the faster convergence 
when the initial condition /(x.r')==/,(x). This comment 
is relevant to the numerical implementation of the PIS 
nicth<xl. and will be discussed below. 

To obtain the PIS appropriate for the dynamic systems 
studied in this paper, it is necessary to be more specific 
on the strocture of the matrix function <2(*)- In particu- 
lar, it will be assumed that the first r rows of Q(,-) are 
zero, that is 



^,^(0=0 for /=1 r;}-\,...Mi ir<n) 



(10) 



and that q,j(*)^0 for al least one j for every i^r+1 n. 

This implies that tlie diffusion matrix G(*) assumes the 
form 



G{')= 



O 




O 
6(-) 



(II) 



j^p{xAx\f)~-^^{m,{x)p{xAx\n\ 



n n 



+ 5 S 2 ^ \g>Ax)pix,t I x\n\. (12) 

^ i^^x ^„i «X(rtXj 

Proceeding in a manner similar to the derivations in 
Ref. [7], it can be shown that the TPD for small values 
of t(-/— r') is given by the following expression, which 
is correct up to terms of order t' 



pixj+T I x',/)-{n5(x,-x/-m,(x')T)} 



.(2iiT)-V|e(x')l4 •cxp{-:^ 2 S 



2t 



f— r+l ^r+1 



{x-x:~mXx')7) [G(x') 'J, ,j-,(Xj-x/-mXx')T)}, 

(13) 

where I G | denotes the determinant of the reduced dif- 
fusion matrix 6, assumed to be positive definite. This 
implies that 1 6 1 > 0. \G~\ denotes the element in po- 
sition ij of the inverse matrix of G. As shown in Ref. 
(71, the expression given by Eq. (13) is not unique, but 
seems to be well suited for our purpose. 

Having obtained an explicit expression for the TPD 
for a .short time step, one Ciui now invoke the Markov 
property. This allows a TPD over a time interval of 
arbitrary length to be expressed in terms of a product of 
short-time TPDs. By dividing a given time interval {t'j) 
into N small time inter\'als of length T^(t—t')/N, it is 
found that (f,=f'+/T, f=f'', x-x*'^'*, r=h. x'=x"") 



(x,^ I x',f')=|-[np(x<'Vj I x"-"V,)dx*"".dr' 



tn II 

(14) 



Similarly, with an initial PDF/(x',f')=H'(x'). the PDF 
f{x,i) will be given by 

/(x,/)=f—fnp(x^\fjx*^".0 ,)H'(x«"Odx*'"."dx'-^ ". 

(15) 

Hence, by combining Eq. (13) withEqs. (14) or (1 5), 
a formal (approximate) solutitin of the FPK equation can 
be written. Equations (14) and (15), which arc often 
referred to as PIS, constitute the core of the numerical 



467 



Volun»e 99. Number 4. July-Augusl 1994 

Journal of Research of the National Institute of Standards and Technology 



solution procedure to be described subsequently. It is 
realized that a numerical solution according to this 
method, automatically provides the evolution in time of 
the (conditional) PDF of the Markov process X(,l) from 
given start conditions in terms of an initial density 
f{x' ,t')=w{x'), including the degenerate case 
/(x',f')-8{x'— Xo), for some starting point Xa. It is also 
worth noting how the PIS relates to the physics of the 
dynamic model, which is expressed through the coeffi- 
cients ffij(») and qiji*), cf. Eq. (3). The evolution in time 
of the PDF as expressed by the PIS, is seen to be directly 
determined by these coefficients in an explicit manner. 
This fact is a very important advantage of the PIS 
method, and reveals its fundamental physical signifi- 
cance. 



3. Numerical Implementation 

In the numerical implementation, the PIS is obtained 
by an iteration process based on the Chapman-Kol- 
mogorov equation expressed as 

p(jc<'Vj I x',n^j'*;jp(x^\tj\x'J-'\tj-t)pix'^-'K 

0_,U',x')dx<'-". (16) 

The discretization of state space for the numerical 
solution makes it appropriate to employ an interpolation 
and smoothing procedure to increase the numerical effi- 
ciency. It was found that application of cubic B-spiines, 
as detaikd in Ref [ 1 7], offered the desired accuracy and 
smoothness for the type of problems considered in this 
paper. This procedure was used as follows. At each time 
step r^|-*r;, p(jr""",/j-i | x',i') is represented as a cubic 
B-spline series in the following manner 



*fi 



p(jc«-",oU'/)-2' 



1^1 



M 



®Bt,(*" "). 



(17) 



where M/=number of grid points for the i'th state vari- 
able x,, {®tiBi,(')}*Jli is a tensor product ba.sis of cubic 
B-splines and {/^ "(Jt|,*",Jta)}JJli is the set of interpola- 
tion coefficients associated with time tj ,. It is assumed 
that each set {Bt.(*)}'^(t J-I,*",n, is a basis of cubic 
B-splincs associated with the knot sequence determined 
by the grid points for the i'th variable Xj. The tensor 
product B-spline is defined by 



®Bt.ix)-tlB,Xx,). 



The representation of p(j:*'"",fj_i \x',t') by B-splincs 
makes it possible to retain high numerical accuracy even 
with a fairly coarse basic grid if p(x''"'V>-i \x',i') is not 
too singular. By substituting from Eq. (17) into Eq. (16), 
Eq. (19) is obtained 

p(xOVj |x',r')=2-2 I^-'\ku"*Jc,) 
\"'\p{x^\tj \x^-'\tj-t) ®5t,(x*^")dx"-". (19) 

J R" J (-1 

It is seen from Eq. (13) that since mj(*) and gy(») are 
not functions of lime t, the TPDs cannot depend on 
absolute time, but only on the time increment. Markov 
processes whose TPDs have this property, are called 
homogeneous. It follows that 

p{x^\tj I x" ",0-,)-p(j:^\t I xf-'\0),y=l ,2 (20) 

which holds for any /;—/;_ i-tsO. 

From Eqs. (19) and (20) it is seen that for a fixed 
value of the time increment t, each of the integrals on 
the right hand side of Eq. (19) need to be calculated only 
once, and can be stored for repeated use. That is, the 
following parameters arc calculated initially and stored 

fii'',-:'r,-|***J/'(*(Vl...,,,.T I :«r^-'>,0) ®B., ('"'")<''"""• 

(21) 

Here, the index /; , i— 1 ,...,«, refers to grid point number 
/( for the state space variable x,. It may be noted here that 
due to the properties of the TPD for small time incre- 
ments T, the tensor fii',:^, has a strongly banded charac- 
ter with the elements decreasing rapidly away from the 
main diagonal /:i=/i,...,jl:„-/„. ThLs has important impli- 
cations for the efficiency of the computer program. Let 
p/f..,^=p(X(VJj^),rj |x',f'). Then Eq. (19) can be rewritten 
as 



1,-1 t.-i 



(22) 



Having calculated the TPD p(x^Vj Ix',/') at the grid 
points by using Eq. (22), a spline interpolation is again 
carried out and a new set of interpolation coefficients 
{P^(k,,...Jc„)}1f^t are calculated. This provides an up- 
dated representation of the TPD for time step /, cf. Eq. 
(17). For each time step, the normalization condition Eq. 
(7) is checked. That is, if 

r r Ml .«« 
J •^•jp(x^'>jj\x\nd^^=l,—J,l^\k A) 



1,-1 VI 



(18) 



fl f B,Xx)dx=qj 

i-l J 



(23) 
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and qjt- 1 .0 within the desired accuracy, then the follow- 
ing replacement is made to restore the correct normal- 
ization. 



Z+Zg{E)+h{Z)~rN{(). 



(28) 



r^Xku'^xr- i-cif mk,,'",k„r 



(24) 



This normalization check and replacement strategy con- 
tributes to producing a ^■ery stable and accurate numer- 
ical procedure. 



4. Examples 

The accuracy and power of the developed PIS prtKe- 
durc wi!! be illustrated by application to specific case 
studies taken from two classes of dynamic models. Both 
models are described by Eq. (3) with n-2 and m-3. 
This implies a two-dimensional state space vector 
X^iXtJCif'^iZJZf. Further, m(-) and e(*) are such that 
'ni{X\JCi)-Xj and ^ij(')-0 fory=l,2,3. Assuming suffi- 
cient restrictions on m{») and Q{'}, cf. Refs. [6,7), X{t) 
becomes a Markov diffusion process. Im'oking Eq. ( 1 3), 
it can be shown that, up to correction terms of order tS 
the associated TPD assumes the form 

P(j:.7 |x',0)=S(jr,-Jc,'-j:2' T>^(jr2,T | x'.O). (25) 

P{x2,t\x',0) is given by the relation 

I 



p{Xi,T\x',0}- 



where 



V2Tr^(jf) 



LAU (-^2-x,'-m,(x')r)^ 
■rfr n 2^(jc')T /• 

(26) 



P(Jc')-2 ivix')' 



(27) 



By combining Eqs. (25) and (26), and applying the 
soluCbn technique described in the previous section, the 
TPD p{x,t \x',i') for large t—t' can be calculated. By 
this, the time evolution of the system when it starts from 
rest, for example, can be studied. The stationary PDF is 
obtained in the limit as t—t'-^°o. For application of the 
PIS method to other problems involving bodi two- and 
three-dimcnsiona! slate space vectors, the reader m^ 
consult Reft. [11-13,18,19]. 

4,1 Example 1 — ^The Caughey Oscillator 

There is a class of dynamic models for which there 
exist an analytical solution for the stationary joint PDF 
of X. A member of this cla.ss m^ be called a Caughey 
oscillator, Ref. [20]. The generic equation of motion for 
this oscillator can be written as 



N{t) denotes a stationary, zero-mean Gaussian white 
noise satisfying E[N(i) A'(r+T)]-5(T), where S(«) de- 
notes Dirac's delta function, /'is a positive constant and 
g{E) is a function of the total energy E=E(Z^) given 
as follows 



E=^Z^+V(Z) 



(29) 



where 



"i 



V(z)- h{s)ds. 



(30) 



For this example "12(2,2) — jgl£(z,j)]-/i(z), and 
we may choose ^2i=^2j=Wi-W2=0. gij-'F and 
dWi(t)-N{t)di. The stationary, joint PDF, denoted by 
Ps{'}, is then determined by the relation, cf, Refs. 
[20,21] 

Ps{z4)-C SKpi- j^ \ g(s)ds}, (31) 



where E'-^z^/2+V{z), and C is a normalization constant 
to ensure a total probability equal to 1.0. 

For the illustration purposes in this paper, we have 
chosen the following special case of Eq. (28) 

Z(i)+2^{t){l+ey^Z\t)+^Z\t)+\xZ*it)]"'}+Zit) 



+AZ^(f)=2\/^(0 



(32) 



with parameters ^, e, and A. 

The stationary PDF only depends on the parameters e 
and A, and the numerical solution for the following set 
of parameter values has been calculated (e,A)=(0.0) 
(Gaussian response), (0, 0.2) (Duffing oscillator) and 
(0.5, 0.1). The calculations were carried out with the 
same number of grid points on both axes in state space, 
aviz., 45. Since the resulting PDFs are actually indepen- 
dent of f , the value ^=0.1 was chosen for the Gaussian 
and Duffing cases, while ^-0.5 was adopted for the last 
case. The time inoements used were 7^.(X)25 s, O.OOl 
s, and 0.02 s, respectively. The total CPU time on a DEC 
station 3 100' was about 5 minutes for each case. In Figs.. 



'Certain commercial etniipment. inslnimcnts, or materials are identi- 
fied in this fopcr 10 spcciiy adequately the experimental procedure. 
Such idenlificalion does not imply recommendalkm or endorsemcM 
by the Nauonal Institute of Standards and Technofogy, nor docs it 
imply that the materials or equipment identified arc necessarily the 
best available for the purpose. 
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1 and 2 are shown the marginal PDFs of the displace- 
ment response for the three case studies considered, 
together with the corresponding analytical solutions. In 
Fig. 3 are giN-en the corresponding analytical and numer- 
ical results for the mean upcrossing rate. It is seen that 
in all three cases the agreement between the numerical 
PIS and the analytical solution is very gixxl over the 
whole range of probability levels given. In fact, the 
accuracy can be retained down to much lower probabil- 
ity lewis (=10" '") at a mtxlcratc increase in computer 
time. 

4^ Example 2 — Parametric and External 
Excitation 

In this example, the response statistics of a nonlinear 
oscillator subjected to both external and parametric ran- 
dom excitation will be illustrated by applying the 
methodology of the paper to two specific case studies. 

The equation of motion of the oscillator is the follow- 
ing 



Z+2a'^+mO]Z+y[Z'+—]Z+c^\\+N2(t)]Z~N,(t). 

(33) 

Here ^, y, and i«Jti are positive constants, Njit ),]"], 2,3, 



arc indepetident Gaussian white noises satisfying 

ElNji,)Ni{i+T)hrf 6{T).>-1.2.3. (34) 

where /} are positive constant.s. For this example it 
is found that m^(^,^)--2f^-■y[^^+^Vt^>,?l-&)^?^, 
qii(zj.)~-2£zrt. qnizj.)--(^zri and qi^izJ.)^!^- 

This model was studied by Dimenttierg [22], who 
showed that when 



Wi?=4f J? 



(35) 



a closed-form cxpressbn for the stationary joint PDF 
can be obtained. It is given as 






where C is a normalization constant and 



n 



.s- 



2^ 



a^n^^aj^^n 2"^ W^ri 



1 



TIF: 



(36) 



(37) 



By this, we have the opportunity to test the accuracy of 
the PIS method for this kind of dynamic model. The 
results of two pitfticular cases will be presented. 




-1 1 

DtiplaccTiKnE 



Fig. 1. IVobabiliiy density function of displacement response for ihe Caughey oscillator in example I. Ana- 
lytical solutions: , e-0. A-0; -•-•, c-0, A-0.2; . e-0.5. A-O.I, Numerical path integral 

solution: +, e-0, A-0; x, £=t), A=0,2; O, »:-0.5, A-fl. I. 
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Fig. 2. Logarithmic plot of the probatiility density function of displacement response for the Cattghey oscillator 
in example 1. Key as in Fig. !. 
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Fig. 3, Mean upcrossing rale of Uisplacemcnt response for the Caughey oscillator in example 1 . Key as in Fig. 
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Case 1: Here the following parameter values were 
used. ^-1, -H).!, Wu-1.0, /7=2.5, /^^O.l, fJ-O.S. 
For the numerical calcu1atk>ns a grid size of 49X49 
points and a time increment 7M).01 s was used. The total 
CPU time on a DEC 3100 work station was 3 min for the 
PIS calculation. The results for the analytical and nu- 
merical solutions are given in Figs. 4-6. In Figs. 4 and 
5 are shown the marginal PDF of the displacement re- 
sponse and in Fig. 6 is shown the correspor>ding mean 
upcrossing rate. 

Case 2: In this case the following set of parameters 
were used. ^.I, yOA, a»o-I.O, /7=5.0, ^-0-2, 
rj-0,3. A grid size of 51X51 points together with a 
time increment 7M>.0l s were chosen. The CPU time 
was about the same as in the previous case. The same 
results as for Case 1 are presented in Figs. 4-6. 



5. Conclusions 

A numerical method for estimating the extreme re- 
sponse of nonlinear oscillators excited by white noise, 
or filtered white noise, has been described. The example 
calculations presented show that the method gives very 
accurate estimates of the required joint PDF. In fact, for 
every example having analytical solution on which the 
method has been tested, complete agreement has been 
found with proper choice of grid size atxl time incre- 
ment in the numerical solution procedure. In the present 
paper, of course, only a few cases can be given. Experi- 
ence with the method indicates that two-dimensional 
problems can be solved routiuely with high accuracy 
requiring a few minutes CPU time on a work statbn 
(DEC station 3100). The solutwn of three-dimensional 
problems requires more care in the sense that computer 
capacity starts to become an issue of importance. In 
such cases the CPU time easily runs into hours. 
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Fig. 4. Probabiliiy density function of displacement response for the oscillaior in example 2, case 1 and 2. 
Analjtical solutions: — . case 1; - ■ - *, case 2. Numerical palh integral solution: +, case I; O. case 2. 
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Fig. 5. Logarithmic plot of the t»-obabiliiy density function of displacement response for the oscillator in 
example 2, case 1 and 2, Key as in Fig, 4. 
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Fig. 6. Mean upcrossing rate of displacemcni response for the oscillalor in example 2. case 1 and 2. Key as in 
Fig. 4. 
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